#pragma rtGlobals=1		// Use modern global access method.

//----------------------------------------------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------------------------------------------

Function steepestDescent(inputMatStr, inputVectStr, startVectStr, iterations)
String inputMatStr, inputVectStr, startVectStr
Variable iterations
Wave inputMat = $inputMatStr, inputVect = $inputVectStr, startVect = $startVectStr

Variable p1 = 0, del, del0, alpha, p2
MatrixOp/O residuals = inputVect -  inputMat x startVect
MatrixOp/O delHolder = residuals^t x residuals 
del = delHolder[0]
del0 = del

DoWindow OptGraph
if(V_flag)

else
		Display/N=OptGraph inputVect, startVect
endif


Do
	MatrixOP/O q_vect = inputMat x residuals
	MatrixOP/O alpha_wave = (residuals^t x q_vect)
	alpha = del/alpha_wave[0]
	//if(del < 6000)
	//	alpha =alpha*10
	//endif
	WaveStats/Q startVect
	
	p2 = 0
	Do
		If(startVect[p2] < (V_sdev - V_adev))
			startVect[p2] = startVect[p2] + .001*alpha*residuals[p2]
		else
			startVect[p2] = startVect[p2] + alpha*residuals[p2]
		endif
		p2+=1
	while(p2 < numpnts(startVect))
	
	p2=0
	Do
		If(startVect[p2] < 0)
			startVect[p2] = 0
		endif
		p2+=1
	while(p2 < numpnts(startVect))
	
	MatrixOp/O residuals = inputVect - (inputMat x startVect)
	MatrixOp/O delHolder = residuals^t x residuals 
	del = delHolder[0]
	if(mod(p1,2) == 0)
	//	Smooth 3,startVect
	endif
	DoUpdate/W=OptGraph
	p1+=1 
while(p1 < iterations && del > .00000000001)

end

//----------------------------------------------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------------------------------------------

Function conjugateGradient(inputMatStr, inputVectStr, startVectStr, iterations)
String inputMatStr, inputVectStr, startVectStr
Variable iterations
Wave inputMat = $inputMatStr, inputVect = $inputVectStr, startVect = $startVectStr

Variable p1 = 0, del, del0, alpha
MatrixOp/O residuals = inputVect -  inputMat x startVect
MatrixOp/O d_vect = residuals
MatrixOp/O delHolder = residuals^t x residuals 
del = delHolder[0]
del0 = del
DoWindow OptGraph
if(V_flag)

else
	Display/N=OptGraph inputVect, startVect
endif

Do
	MatrixOP/O q_vect = inputMat x d_vect
	MatrixOP/O alpha_wave = (d_vect^t x q_vect)
	alpha = del/alpha_wave[0]
	startVect = startVect + alpha*d_vect
	Variable p2 = 0
	
	Do
		If(startVect[p2] < 0)
			startVect[p2] = 0
		endif
		p2+=1
	while(p2 < numpnts(startVect))
	
	MatrixOp/O residuals = inputVect - (inputMat x startVect)
	Variable delOld = del
	MatrixOp/O delHolder = residuals^t x residuals 
	del = delHolder[0]
	Variable betaValue = del/delOld
	d_vect = residuals + betaValue*d_vect
	
	if(mod(p1,5) == 0)
		//Smooth 2,startVect
	endif
	
	DoUpdate/W=OptGraph
	p1+=1 
while(p1 < iterations && del > .00000000001)

end

//----------------------------------------------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------------------------------------------

Function steepestMaxEnt(inputMatStr, inputVectStr, startVectStr, iterations, lambda)
String inputMatStr, inputVectStr, startVectStr
Variable iterations,lambda
Wave inputMat = $inputMatStr, inputVect = $inputVectStr, startVect = $startVectStr

Variable p1 = 0, del, del0, alpha, var, p2
MatrixOP/O vararr = powR((inputVect - mean(inputVect)),2)/(numPoints(inputVect)-1)
var = sum(vararr)
MatrixOp/O d2chi = inputMat x (inputVect - inputMat x startVect)

p2=0
Do
	If(startVect[p2] <= 0)
		startVect[p2] = 10^-5
	endif
	p2+=1
while(p2 < numpnts(startVect))

Make/O/N=(numpnts(startVect)) oneWave = 1

Variable coefficent = (2*lambda / var)
d2chi = coefficent * d2chi

MatrixOp/O gradQ = oneWave + ln(startVect/(sum(startVect))) +  d2chi
MatrixOp/O alphaWave = (startVect^t x gradQ)/(startVect^t x startVect)
alpha = alphaWave[0]

DoWindow OptGraph
if(V_flag)

else
		Display/N=OptGraph inputVect, startVect
endif


Do
	MatrixOp/O gradQ_old = gradQ
	MatrixOp/O startVect_old = startVect
	
	startVect = startVect - 1/alpha*gradQ
	MatrixOp/O d2chi = inputMat x (inputVect - inputMat x startVect)
	d2chi = coefficent * d2chi
	p2=0
	Do
		If(startVect[p2] <= 0)
			startVect[p2] = 10^-5
		endif
		p2+=1
	while(p2 < numpnts(startVect))
	MatrixOp/O gradQ = (oneWave+ ln(startVect/sum(startVect)) + d2chi)
	MatrixOp/O alphaWave = ((startVect-startVect_old)^t x (gradQ-gradQ_old))/((startVect-startVect_old)^t x (startVect-startVect_old))
	alpha = alphaWave[0]
	

	DoUpdate/W=OptGraph
	p1+=1
while(p1 < iterations)

end

//----------------------------------------------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------------------------------------------

Function gradientMaxEnt(inputMatStr, inputVectStr, startVectStr, iterations, lambda)
String inputMatStr, inputVectStr, startVectStr
Variable iterations,lambda
Wave inputMat = $inputMatStr, inputVect = $inputVectStr, startVect = $startVectStr

Variable p1 = 0, del, del0, alpha, var, p2
MatrixOp/O residuals = inputVect -  inputMat x startVect
MatrixOp/O d_vect = residuals
MatrixOp/O delHolder = residuals^t x residuals 
del = delHolder[0]
del0 = del
MatrixOP/O vararr = powR((startVect - mean(startVect)),2)/(numPoints(startVect)-1)
var = sum(vararr)
MatrixOp/O d2chi = inputMat^t x (inputVect - inputMat x startVect)

p2=0
Do
	If(startVect[p2] < 10^-5)
		startVect[p2] = 10^-5
	endif
	p2+=1
while(p2 < numpnts(startVect))

Make/O/N=(numpnts(startVect)) oneWave = 1

Variable coefficent = (2*lambda / var)
d2chi = coefficent * d2chi
MatrixOp/O gradQ = oneWave + ln(startVect/(sum(startVect))) +  d2chi
MatrixOp/O dir1 = -gradQ
MatrixOp/O alphaWave =(dir1^t x (inputMat x dir1))
alpha = del/alphaWave[0]


DoWindow OptGraph
if(V_flag)

else
		Display/N=OptGraph inputVect, startVect
endif

	Execute "progresspanel()"


Do
	ValDisplay valdisp0, win=ProgressPanel,  limits={0,iterations,0}, barmisc={0,0}, value=_NUM:p1
	DoUpDate /W=ProgressPanel /E=1
	
	MatrixOp/O gradQ_old = gradQ
	MatrixOp/O startVect_old = startVect
	Variable del_old = del
	
	startVect = startVect - alpha*dir1
	MatrixOp/O residuals = inputVect -  inputMat x startVect
	MatrixOp/O delHolder = residuals^t x residuals 
	del = delHolder[0]
	Variable betta = del/del_old
	MatrixOp/O d2chi = inputMat^t x (inputVect - inputMat x startVect)
	MatrixOP/O vararr = powR((startVect - mean(startVect)),2)/(numPoints(startVect)-1)
	var = sum(vararr)
	coefficent = (2*lambda / var)
	d2chi = coefficent * d2chi
	p2=0
	Do
		If(startVect[p2] < 10^-5)
			startVect[p2] = 10^-5
		endif
		p2+=1
	while(p2 < numpnts(startVect))
	MatrixOp/O gradQ = (oneWave+ ln(startVect/sum(startVect)) + d2chi)
	MatrixOp/O alphaWave =(dir1^t x (inputMat x dir1))
	alpha = del/alphaWave[0]

	MatrixOp/O dir1 = -gradQ + betta*dir1
	DoUpdate/W=OptGraph
	p1+=1
while(p1 < iterations)
KillWindow ProgressPanel
end

//----------------------------------------------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------------------------------------------

